/**
 * Copyright (C) 2025  粤港澳大湾区（广东）国创中心
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along
 *  with this program; if not, see <https://www.gnu.org/licenses/>.
 */
#include "NctiGmshPubfun.h"
#include "GmshDefines.h"

#include "NctiDisplayModelBody.h"

#include "NctiMeshManager.h"
#include "NctiOCCBody.h"

#include "gmsh.h"
#include <NctiTime.h>
#include <map>
#include <NctiOcTree.h>
#ifdef _MSC_VER
#define NOMINMAX
#define _CRT_SECURE_NO_WARNINGS
#include <windows.h>
#endif // !_MSC_VER	
#include "NctiDebug.h"
#include <algorithm>

#include "NctiGEFunction.h"
#include "NctiDisplayDocCustomContainer.h"
#include "NctiGeometryFunction.h"
#include "NctiBaseUtility.h"
#include "NctiString.h"
#include "NctiBaseModel.h"
#include "NctiDisplayObjectLifeCycleManager.h"
#include "NctiBaseView.h"

using namespace NCTI;
using namespace gmsh;

static bool convert(int gmeshType, NctiElementType&ret) {

    ret = NctiElementType (- 1);
    switch (gmeshType)
    {
    case MSH_PNT:
        ret = NctiElementType::Point1;
        break;
    case MSH_LIN_2:
        ret = NctiElementType::Edge2;
        break;
    case MSH_LIN_3:
        ret = NctiElementType::Edge3;
        break;
    case MSH_TRI_3:
        ret = NctiElementType::Triangle3;
        break;
    case MSH_TRI_6:
        ret = NctiElementType::Triangle6;
        break;
    case MSH_QUA_4:
        ret = NctiElementType::Quadrangle4;
        break;
    case MSH_QUA_8:
        ret = NctiElementType::Quadrangle8;
        break;
    case MSH_QUA_9:
        ret = NctiElementType::Quadrangle9;
        break;
    case MSH_TET_4:
        ret = NctiElementType::Tetrahedron4;
        break;
    case MSH_TET_10:
        ret = NctiElementType::Tetrahedron10;
        break;
    case MSH_HEX_8:
        ret = NctiElementType::Hexahedron8;
        break;
    case MSH_HEX_20:
        ret = NctiElementType::Hexahedron20;
        break;
    case MSH_HEX_27:
        ret = NctiElementType::Hexahedron27;
        break;
    case MSH_PRI_6:
        ret = NctiElementType::Prism6;
        break;
    case MSH_PRI_18:
        ret = NctiElementType::Prism18;
        break;
    case MSH_PYR_5:
        ret = NctiElementType::Pyramid5;
        break;
    case MSH_PYR_13:
        ret = NctiElementType::Pyramid13;
        break;
    case MSH_PYR_14:
        ret = NctiElementType::Pyramid14;
        break;
    default:
        break;
    }
    return ret>=0;
}

NctiGmeshHelper::NctiGmeshHelper(Ncti_Mesh_Option op)
{
    m_op = op;
    gmsh::initialize();
    model::add("test");
}

NctiGmeshHelper::~NctiGmeshHelper()
{
    model::remove();
    gmsh::finalize();
}

int RunGMsh(bool volume)
{
	int rt = -1;
#ifdef _MSC_VER
	__try
	{
#endif // !_MSC_VER	
		if (volume)
		{
			model::mesh::generate(3);
		}
		else
		{
			model::mesh::generate(2);
		}
		//std::cout << "result:" << rt << std::endl;
        rt = 1;

#ifdef _MSC_VER	
	}
	__except (EXCEPTION_EXECUTE_HANDLER)
	{
        std::cout << "GMsh Error!" << rt << std::endl;
		rt = 0;
	}
#endif // !_MSC_VER	
    return rt;
}

bool NctiGmeshHelper::CalculaeSurface()
{
    //option::setNumber("General.Terminal", 1);
    option::setNumber("Mesh.MeshSizeMax", m_op.m_size);
    option::setNumber("Mesh.MeshSizeMin", m_op.m_size);
    option::setNumber("Mesh.AlgorithmSwitchOnFailure", 1);
    option::setNumber("Mesh.MeshSizeFromCurvature", m_op.m_autoCurvature);
    if (m_op.m_angleTol > 0)
    {
        option::setNumber("Mesh.AngleToleranceFacetOverlap", m_op.m_angleTol);
    }
    //option::setNumber("Mesh.MinimumCirclePoints", 6);
    //option::setNumber("Mesh.MeshSizeMin", op.m_size * 1.1);
    
    if (m_op.m_high)
        option::setNumber("Mesh.ElementOrder", 2);
    //option::setNumber("Mesh.MeshSizeFromCurvature", op.m_auto);

	if (m_op.m_bParallel)
	{
		int numThreads = std::thread::hardware_concurrency() - 1;
		// 设置并行参数
		gmsh::option::setNumber("Mesh.MaxNumThreads1D", numThreads);  // 1D网格线程数
		gmsh::option::setNumber("Mesh.MaxNumThreads2D", numThreads);  // 2D网格线程数
		//gmsh::option::setNumber("Mesh.MaxNumThreads3D", numThreads);  // 3D网格线程数
	}
    try
    {
        NctiTime time("Gmesh generate 2d", true);
        if (RunGMsh(false) != 1)
            return false;
        //model::mesh::generate(2);
    }
    catch (...)
    {
        std::string str;
        logger::getLastError(str);
        std::cout << str;
        return false;
    }
    return true;
}

bool NctiGmeshHelper::CalculateVolume()
{
    option::setNumber("Mesh.MeshSizeMax", m_op.m_size);
    option::setNumber("Mesh.MeshSizeMin", m_op.m_size);
    option::setNumber("Mesh.AlgorithmSwitchOnFailure", 1);
    option::setNumber("Mesh.MeshSizeFromCurvature", m_op.m_autoCurvature);
    if (m_op.m_angleTol>0)
    {
        option::setNumber("Mesh.AngleToleranceFacetOverlap", m_op.m_angleTol);
    }
    //if (m_op.m_qualityOpt>0)
    {
        option::setNumber("Mesh.Optimize", 1);
        //option::setNumber("Mesh.MinQuality", m_op.m_qualityOpt);
    }
    
    //option::setNumber("Mesh.MinimumCirclePoints", 6);
    //option::setNumber("Mesh.MeshSizeMin", op.m_size * 1.1);

    if (m_op.m_high)
        option::setNumber("Mesh.ElementOrder", 2);

    //option::setNumber("Geometry.OCCSewFaces", 1);
    //option::setNumber("Geometry.Tolerance", 1e-4);
	//model::occ::healShapes(dimTagsEx);
    if (m_op.m_bParallel)
    {
        int numThreads = std::thread::hardware_concurrency() - 1;
        // 设置并行参数
        gmsh::option::setNumber("Mesh.MaxNumThreads1D", numThreads);  // 1D网格线程数
        gmsh::option::setNumber("Mesh.MaxNumThreads2D", numThreads);  // 2D网格线程数
        gmsh::option::setNumber("Mesh.MaxNumThreads3D", numThreads);  // 3D网格线程数
    }
	try
	{
		NctiTime time("Gmesh generate 3d", true);
        if (RunGMsh(true) != 1)
            return false;
        //model::mesh::generate(3);
    }
    catch (...)
    {
        std::string str;
        logger::getLastError(str);
        std::cout << str;
        return false;
    }
    return true;
}

void setNonQuadMeshSize(int surfaceTag, double targetMeshSize) {
	try {
		// 1. 获取表面的所有顶点（用于设置网格尺寸的参考点）
		std::vector<std::pair<int, int>> surfEntities = { {2, surfaceTag} };
		std::vector<std::pair<int, int>> vertices;
		gmsh::model::getBoundary(surfEntities, vertices, false, false, true);  // 递归到顶点（dim=0）

		//// 2. 提取顶点标签，设置每个顶点的网格尺寸（控制三角形边长）
		//std::vector<int> vertexTags;
		//for (const auto& [dim, tag] : vertices) {
		//	if (dim == 0) {  // 仅保留顶点（dim=0）
		//		vertexTags.push_back(tag);
		//	}
		//}
		// 3. 设置顶点的网格尺寸（Gmsh 会自动插值生成均匀的三角形网格）
		gmsh::model::mesh::setSize(vertices, targetMeshSize);
		//std::cout << "非四边形 Surface " << surfaceTag << " 网格尺寸设置为：" << targetMeshSize << std::endl;
	}
	catch (const std::exception& e) {
		std::cerr << "设置非四边形 Surface " << surfaceTag << " 网格尺寸失败：" << e.what() << std::endl;
	}
}

void FullQuad(NctiDisplayDocCustomContainer* i_pCustomContainer, NctiOCCBody* pOCCBody, const std::map<std::pair<int, int>, int>& mapParts, const Ncti_Double& meshSize)
{
	//gmsh::option::setNumber("Mesh.Algorithm", 5);
	std::vector<std::pair<int, int>> surfaces;
	gmsh::model::getEntities(surfaces, 2);
	std::map<int, int> mapEdgeNumPt;

	NctiBaseView* pNctiBaseView = i_pCustomContainer->GetLCM()->GetBaseModel()->GetActiveBaseView();
	NctiRenderSegmentItem* pSubSegment = pNctiBaseView->GetPreviewKey("FullQuadtext");
    pSubSegment->GetVisibilityControl()->SetEverything(true);
    pSubSegment->Flush();
	for (const auto& surface : surfaces)
	{
		int dim = surface.first;
		int tag = surface.second;
		try
		{
            gmsh::vectorpair loops;
			gmsh::model::getBoundary({ surface/*{2, tag}*/ }, loops, true, false, false);
			std::sort(loops.begin(), loops.end());
            loops.erase(std::unique(loops.begin(), loops.end()), loops.end());
			std::vector<int> boundaryCurves;
			std::vector<int> nctiTags;
            for (int i = 0; i < loops.size(); i++)
			{
				if (loops[i].first != 1)//只保留曲线
                    continue;
				auto it = mapParts.find(loops[i]);
				if (mapParts.end() == it)
					continue;
				boundaryCurves.push_back(std::abs(loops[i].second));
                nctiTags.push_back(it->second);
            }
            
			// 1. 获取当前面的所有边界曲线（修正后的数据类型，确保获取正确的曲线标签）
			//std::vector<int> boundaryCurves = getSurfaceBoundaryCurves(tag);
			int numCurves = boundaryCurves.size();

			// 仅对四边形表面（4条边）使用 Transfinite 结构化网格
			if (numCurves != 4) {
				// 步骤1：设置均匀的三角形网格尺寸（关键：确保重组后四边形均匀）
				setNonQuadMeshSize(tag, meshSize);

				// 步骤2：启用重组（三角形→四边形）
				gmsh::model::mesh::setRecombine(dim, tag);

				// 步骤3：可选优化：设置重组参数（提高网格质量，4.11.1 兼容）
				// - 最小重组角度（默认 10°，可调整为 15° 减少畸形）
				//gmsh::option::setNumber("Mesh.Recombine.MinAngle", 30.0);
				// - 允许非凸多边形重组（默认关闭，开启后支持凹多边形）
				gmsh::option::setNumber("Mesh.Recombine.AllowNonConforming", 1);

				//std::cout << "非四边形 Surface " << tag << " → 三角形重组四边形网格" << std::endl;
				//std::cout << "Surface " << tag << " 不是四边形（含 " << numCurves << " 条边），跳过 Transfinite，仅重组网格" << std::endl;
				//gmsh::model::mesh::setRecombine(dim, tag);  // 三角形重组为四边形
				continue;
			}
            else
			{
				int curve1 = boundaryCurves[0];  // 边1（如底边）
				int curve2 = boundaryCurves[1];  // 边2（如右边）
				int curve3 = boundaryCurves[2];  // 边3（如顶边，与边1成对）
				int curve4 = boundaryCurves[3];  // 边4（如左边，与边2成对）

                int tag1 = nctiTags[0];
                int tag2 = nctiTags[1];
                int tag3 = nctiTags[2];
                int tag4 = nctiTags[3];

				NctiPosition pt1_S, pt1_E, pt2_S, pt2_E, pt3_S, pt3_E, pt4_S, pt4_E;
				NctiVector vecS, vecE;
				NctiGeometryFunction::get_edge_end_points(i_pCustomContainer,
                    pOCCBody, &tag1, pt1_S, pt1_E, vecS, vecE);
				NctiGeometryFunction::get_edge_end_points(i_pCustomContainer,
                    pOCCBody, &tag2, pt2_S, pt2_E, vecS, vecE);
				NctiGeometryFunction::get_edge_end_points(i_pCustomContainer,
                    pOCCBody, &tag3, pt3_S, pt3_E, vecS, vecE);
				NctiGeometryFunction::get_edge_end_points(i_pCustomContainer,
                    pOCCBody, &tag4, pt4_S, pt4_E, vecS, vecE);

                double d1 = DBL_MAX, d2 = DBL_MAX, d3 = DBL_MAX, d4 = DBL_MAX;
                d1 = distance_to_point(pt1_S, pt3_S);
                d2 = distance_to_point(pt1_S, pt3_E);
                d3 = distance_to_point(pt1_E, pt3_S);
                d4 = distance_to_point(pt1_E, pt3_S);
                if (d1 < 1e-3 || d2 < 1e-3 || d3 < 1e-3 || d4 < 1e-3)
				{
					d1 = distance_to_point(pt1_S, pt2_S);
					d2 = distance_to_point(pt1_S, pt2_E);
					d3 = distance_to_point(pt1_E, pt2_S);
					d4 = distance_to_point(pt1_E, pt2_S);

                    if (d1 < 1e-3 || d2 < 1e-3 || d3 < 1e-3 || d4 < 1e-3)
					{
						int val = curve4; curve4 = curve3; curve3 = val;
						val = tag4; tag4 = tag3; tag3 = val;
                    }
                    else
                    {
						int val = curve2; curve2 = curve3; curve3 = val;
						val = tag2; tag2 = tag3; tag3 = val;
                    }
                }

                int numDivisionsY, numDivisionsX;
                auto iteg = mapEdgeNumPt.find(tag1);
                if (iteg == mapEdgeNumPt.end())
                {
                    iteg = mapEdgeNumPt.find(tag3);
                    if (iteg == mapEdgeNumPt.end())
                    {
                        double dlen1 = 0;
                        NctiGEFunction::get_edge_length(i_pCustomContainer, pOCCBody, &tag1, dlen1);

                        double dlen3 = 0;
                        NctiGEFunction::get_edge_length(i_pCustomContainer, pOCCBody, &tag3, dlen3);
                        int numDiv1 = dlen1 / meshSize + 1;
                        int numDiv3 = dlen3 / meshSize + 1;
                        numDivisionsX = std::max(numDiv1, numDiv3) + 1;
                        mapEdgeNumPt[tag1] = numDivisionsX;
                        mapEdgeNumPt[tag3] = numDivisionsX;

						Ncti_Char ch[256];
						NctiString::LocaleSprintf(ch, "%s",  NctiString::int_to_string(numDivisionsX).c_str());
						std::string str(ch);

                        //
                        NctiPosition midPt;
						NctiGeometryFunction::get_edge_mid_point(i_pCustomContainer, pOCCBody, &tag1, midPt);
                        pSubSegment->InsertText(str.c_str(), midPt);
						NctiGeometryFunction::get_edge_mid_point(i_pCustomContainer, pOCCBody, &tag3, midPt);
                        pSubSegment->InsertText(str.c_str(), midPt);
                    }
                    else
                    {
                        numDivisionsX = iteg->second;
                        mapEdgeNumPt[tag1] = numDivisionsX;
                    }
                }
                else
				{
					numDivisionsX = iteg->second;
					mapEdgeNumPt[tag3] = numDivisionsX;
                }
				iteg = mapEdgeNumPt.find(tag2);
				if (iteg == mapEdgeNumPt.end())
				{
					iteg = mapEdgeNumPt.find(tag4);
                    if (iteg == mapEdgeNumPt.end())
                    {
                        double dlen2 = 0;
                        NctiGEFunction::get_edge_length(i_pCustomContainer, pOCCBody, &tag2, dlen2);

                        double dlen4 = 0;
                        NctiGEFunction::get_edge_length(i_pCustomContainer, pOCCBody, &tag4, dlen4);
                        int numDiv2 = dlen2 / meshSize + 1;
                        int numDiv4 = dlen4 / meshSize + 1;
						numDivisionsY = std::max(numDiv2, numDiv4) + 1;
						mapEdgeNumPt[tag2] = numDivisionsY;
						mapEdgeNumPt[tag4] = numDivisionsY;

						Ncti_Char ch[256];
						NctiString::LocaleSprintf(ch, "%s", NctiString::int_to_string(numDivisionsY).c_str());
						std::string str(ch);
						//
						NctiPosition midPt;
						NctiGeometryFunction::get_edge_mid_point(i_pCustomContainer, pOCCBody, &tag2, midPt);
                        pSubSegment->InsertText(str.c_str(), midPt);
						NctiGeometryFunction::get_edge_mid_point(i_pCustomContainer, pOCCBody, &tag4, midPt);
                        pSubSegment->InsertText(str.c_str(), midPt);
                    }
                    else
                    {
                        numDivisionsX = iteg->second;
                        mapEdgeNumPt[tag2] = numDivisionsX;
                    }
				}
                else
				{
					numDivisionsX = iteg->second;
					mapEdgeNumPt[tag4] = numDivisionsX;
                }

                //numDivisionsX = 5;
                //numDivisionsY = 5;
                std::cout << numDivisionsX << "   " << numDivisionsY << std::endl;

				setNonQuadMeshSize(tag, meshSize);
				// 4. 对边设置相同点数（均匀分布），确保 divisions 一致
				gmsh::model::mesh::setTransfiniteCurve(curve1, numDivisionsX);
				gmsh::model::mesh::setTransfiniteCurve(curve3, numDivisionsX);  // 对边点数相同
				gmsh::model::mesh::setTransfiniteCurve(curve2, numDivisionsY);
				gmsh::model::mesh::setTransfiniteCurve(curve4, numDivisionsY);  // 对边点数相同

				// 设置Transfinite曲面
				gmsh::model::mesh::setTransfiniteSurface(tag);

				// 设置Recombine
				gmsh::model::mesh::setRecombine(dim, tag);
            }
		}
		catch (const std::exception& e) {
			std::cerr << "Error processing surface " << tag << ": " << e.what() << std::endl;
		}
	}
}

bool NctiGmeshHelper::ImportModelOCC(NctiDisplayDocCustomContainer* i_pCustomContainer, NctiOCCBody* pOCCBody)
{
    if (!pOCCBody|| !pOCCBody->GetShape())
    {
        return false;
    }

    std::vector<std::pair<int, int>> dimTagsEx;
    auto shape = pOCCBody->GetShape();
    model::occ::importShapesNativePointer(shape, dimTagsEx, true);
    dimTagsEx.clear();
    model::occ::getEntities(dimTagsEx);

    for (int i=0;i<dimTagsEx.size();++i)
    {
        auto& item = dimTagsEx[i];
        int g = -1;
        pOCCBody->get_persistent_from_dimtag(g, item);
        m_parts[item] = g;
    }
    //model::occ::healShapes();  //自动修复几何
    model::occ::synchronize();
    
    if (NctiDebug::getOutputState())
    {
        gmsh::option::setNumber("General.Verbosity", 1);//0=只显示错误，1=错误和警告，2=信息，3=详细信息  打印很耗时！！！！
    } 
    else
    {
        gmsh::option::setNumber("General.Terminal", 0);   // 禁用终端输出
	}
	if (m_op.m_type != 0)
	{
		if (m_op.m_type == 2)
		{
			// 2D算法：Quasi-structured Quad
			gmsh::option::setNumber("Mesh.Algorithm", 11);
			//gmsh::option::setNumber("Mesh.Algorithm", 6);
			// 3D算法：Delaunay
			gmsh::option::setNumber("Mesh.Algorithm3D", 1);
			// 2D重组算法：Simple Full-Quad
			gmsh::option::setNumber("Mesh.RecombinationAlgorithm", 2);
			//“Recombine all triangular meshes”
			gmsh::option::setNumber("Mesh.RecombineAll", 1);
			// 细分算法：None
			gmsh::option::setNumber("Mesh.SubdivisionAlgorithm", 0);
			// 单元阶数：1
			//gmsh::option::setNumber("Mesh.ElementOrder", 1);

            ////Frontal-Delaunay算法需要独立设置边的点数
			//gmsh::option::setNumber("Mesh.Algorithm", 6);
            //FullQuad(i_pCustomContainer, pOCCBody, m_parts, m_op.m_size);
		}
		else
		{
			option::setNumber("Mesh.RecombineAll", 1);
			//0: simple, 1: blossom, 2: simple full-quad, 3: blossom full-quad
			option::setNumber("Mesh.RecombinationAlgorithm", 0);
			//if(op.m_type == 2){//all quad
			//    option::setNumber("Mesh.MeshSizeMax", op.m_size*2);
			//    option::setNumber("Mesh.MeshSizeMin", op.m_size*2);
			//    //0: none, 1 : all quadrangles, 2 : all hexahedra, 3 : barycentric
			//    option::setNumber("Mesh.SubdivisionAlgorithm", 1);
			//} 
		}
	}
    else
	{
		//gmsh::option::setNumber("Mesh.OptimizeNetgen", 1);        
    }
    return true;
}

bool NctiGmeshHelper::CopyResult(
    std::vector<Ncti_Float>& nodePositions,
    std::vector<NctiElementType>& elementType,
    std::vector<Ncti_UInteger>& elementTag,
    std::vector<Ncti_UInteger>& elementData)
{
    NctiLinearMap map;
    std::vector<Ncti_Size> nodeTags;
    std::vector<Ncti_Double> nodeCoords;
    std::vector<Ncti_Double> nodeParams;
    gmsh::model::mesh::getNodes(nodeTags, nodeCoords, nodeParams);

    nodePositions.clear();
    elementType.clear();
    elementData.clear();
    elementTag.clear();
 
    for (int i = 0; i < nodeTags.size(); ++i)
    {
        map[nodeTags[i]] = i;
        auto pos = nodeCoords.data() + 3 * i;
        nodePositions.push_back(pos[0]);
        nodePositions.push_back(pos[1]);
        nodePositions.push_back(pos[2]);
    }

    if (nodeTags.size() > 0)
    {
        std::vector<int> elemTypes;
        std::vector<std::vector<std::size_t> > elemTags;
        std::vector<std::vector<std::size_t> > elemNodeTags;
        for (auto& [dt, tag] : m_parts)
        {
            gmsh::model::mesh::getElements(
                elemTypes, elemTags, elemNodeTags, dt.first, dt.second);
            for (int i = 0; i < elemTypes.size(); ++i)
            {
                NctiElementType tp;
                if (convert(elemTypes[i], tp))
                {
                    auto& nodes = elemNodeTags[i];
                    int count = elemTags[i].size();
                    int numberNodes = nodes.size() / count;

                    for (int j = 0; j < count; ++j) {
                        auto data = nodes.data() + j * numberNodes;
                        elementType.push_back(tp);
                        elementTag.push_back(tag);
                        for (int k = 0; k < numberNodes; ++k)
                            elementData.push_back(map[data[k]]);
                    }
                }
            }
        }
    }
    return nodePositions.size()>0&&elementTag.size()>0;
}

